Gaussian operator bases for correlated fermions 



in 
o 
o 

(N 
> 

o 



> 
O 

o 



in ■ 
o : 

Oh: 

a ■ 
^3 : 



X 



J. F. Corney and P. D. Drummond 
ARC Centre of Excellence for Quantum-Atom Optics, 
University of Queensland, Brisbane 4072, Queensland, Australia. 
(Dated: 1st February 2008) 

We formulate a general multi-mode Gaussian operator basis for fermions, to enable a positive 
phase-space representation of correlated Fermi states. The Gaussian basis extends existing bosonic 
phase-space methods to Fermi systems and thus enables first-principles dynamical or equilibrium 
calculations in quantum many-body Fermi systems. We prove the completeness and positivity of 
the basis, and derive differential forms for products with one- and two-body operators. Because 
the basis satisfies fermionic superselection rules, the resulting phase space involves only c-numbers, 
without requiring anti-commuting Grassmann variables. 



I. INTRODUCTION 

In this paper we address the issue of how to represent 
highly correlated fermionic states, for the purposes of ef- 
ficient calculations in fermionic many-body physics. To 
this end, we introduce a normally ordered Gaussian op- 
erator basis for fermionic density operators. With this 
basis, earlier phase-space techniques used to represent 
atomic transitions 1, 2] can be extended to general Fermi 
systems. 

The analogous phase-space methods for bosons were 
first introduced over a classical phase spacejj ^ IE l^i 
which was subsequently extended to a non-classical phase 
space 0, 1^ and then to a complete Gaussian basis 0. 
Just as with bosons, the Gaussian basis enables a rep- 
resentation of arbitrary fermionic density operators as a 
positive distribution over a generalised phase space. 

We investigate both number-conserving basis sets, 
which generalize the usual thermal states used to de- 
scribe Fermi gases, and non-number-conserving basis 
sets, which generalize the BCS states found in supercon- 
ductivity theory^3- The latter case comprises the most 
general Gaussian basis. 

We concentrate on the foundational issues of the Gaus- 
sian representation method, proving three central results: 

• the basis is complete, 

• the distribution can always be chosen positive, 

• all two-body operators map to a second-order dif- 
ferential form. 

From these results, it follows that positive-definite 
Fokker-Planck equations ^3 exist for many-body 
fermionic systems. Such Fokker-Planck equations enable 
first-principles stochastic numerical simulation methods, 
either in real time or at finite temperature. They can 
also be used to obtain novel types of perturbation theory 
using stochastic diagram techniques l^. 

The formulation of the resultant phase-space simula- 
tion methods, which, for example, can be ap plie d to the 
Hubbard modeljQ, will be given elsewhere jl5i|. How- 
ever, we note these methods have similarities to the aux- 
iliary field method [l3l| used to simplify fermionic path 



integrals, except that here the Gaussian basis is used to 
expand the fermionic states themselves, rather than a 
path integral, which has advantages in terms of giving a 
greater physical understanding and fewer restrictions in 
the resulting applications. 

To begin, we establish in Sec. the definition of a 
Gaussian operator in unnormalised form, and follow this 
in Sec. 11111 with elementary examples of one- and two- 
mode Gaussians in order to illustrate the basic structure 
and the physics underlying it. The two- mode examples 
include the density operators for thermal and squeezed 
fermionic states. In Sec. IIVI we introduce some nota- 
tion that is convenient for discussing general Gaussian 
operators in a multimode system. 

The main part of the paper (Section^ introduces the 
most general class of normalized Gaussian operators and 
discusses those of its properties that are relevant to its 
use as a basis set. This section also gives the proofs of the 
completeness and positivity properties, which are estab- 
lished by expansion in number-state projectors. Other 
properties are proved by use of Fermi coherent states 
(Appendix ^) , and these latter proofs are given in Ap- 
pendix 



II. DEFINITIONS 

To define Gaussian operators for a given fermionic sys- 
tem, we first decompose it into a set of M orthogonal 
single-particle modes, or orbitals. With each of these 
modes, we associate creation and annihilation operators 
bl and bj, with anticommutation relations 



bj,bk 



= 



(2.1) 



where j,k — I...M. Thus, b is a column vector of the 
M annihilation operators, and is a row vector of the 
corresponding creation operators. 

We define a Gaussian operator to be a normally or- 
dered, Gaussian form of annihilation and creation oper- 
ators. Like a complex number Gaussian, the operator 
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Gaussian is an exponential of a quadratic form, with the 
exponential defined by its series representation. We fol- 
low the standard convention that fermionic normal order- 
ing includes a sign-change for every operator exchange, 
so that '■ bkbj := —bjbk- Normal ordering is utilized here 
because it allows us to most directly use Fermi coherent 
state methodsflii involving Grassmann algebrap^[l8l|. 
which are described in greater detail in the Appendices. 
We note that from particle-holc symmetry, it is also pos- 
sible to obtain an entirely analogous anti-normally or- 
dered representation, which simply exchanges the roles 
of particles and holes in the formalism. 

The most general Gaussian form is a cumbersome ob- 
ject to manipulate, unless products of odd numbers of 
operators are excluded. Fortunately, restricting the set 
of Gaussians to those containing only even products can 
be physically justified on the basis of superselection rules 
for fermions. Because it is constructed from pairs of oper- 
ators, this type of Gaussian operator contains no Grass- 
mann variables[l^[l3. 

In order to relate these quadratic expressions to the 
usual Dirac state-vector notation in the examples that 
follow, we define number states as: 



\ni) <Si \n2) 



= b 
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M 



|0)m , (2-2) 



where n = (ni, ^2, . • . , f^M) is a vector of integer occu- 
pations and where |0)^^ is the M-mode fermion vacuum 
state. 

We also recall some elementary identities which relate 
quadratic forms in fermion operators to state projection 
operators. For individual operators, one has the well- 
known identities: 



bj = 10). (1 



3 \ Ij 



ID 



' 3 ^ '3 



(2.3) 



Hence, the elementary identities for quadratic products 
are as follows, for i < j : 



b^b^ - 


ii>,:(ii,ni^- 
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bib, = 












bjbi = 


|Oz,0,)(l„l,| 





(2.4) 



where Ifc is the identity operator for the fc— th mode. 

Expanding quantum states in terms of an underlying 
basis set is a widespread procedure in quantum mechan- 
ics. However, using the general Gaussian operators in 
this capacity is nonstandard in several respects. First, 
one is expanding a quantum density operator over a 



basis set that includes mixed as well as pure quantum 
states. Second, the Gaussian operator basis states are 
not orthogonal. Third, we also allow the basis set to 
include operators which are not themselves density oper- 
ators. These additional degrees of freedom prove useful 
in obtaining the requisite completeness properties. We 
now give the mathematical form of these basis operators, 
in unnormalised form; normalized forms Aj\/, for which 
Tr ( Aj\/) = 1, will be introduced in SecM 



A. Number-conserving Gaussians 

We first consider the special case obtained when the 
quadratic is of form blbj, so that it is number-conserving. 
The most general, number-conserving Gaussian operator 
can be written as: 



: exp 

M 

: Y[ exp -b\lJ.tjb] 

M 

: n [^'b\^^^3^3 



(2.5) 



where /x is a complex M x M matrix. The last result 
follows because normally ordered products in which the 
same operator appears more than once are zero, from the 
anti-commutation relations in Eq. (|2.1|l . For the special 
case that /x is a hermitian matrix, the Gaussian operator 
is just the density operator of a thermal state. If the 
eigenvalues of fi are also either or 1, then the Gaussian 
corresponds to a Slater determinant (i.e. a product of 
single-mode number states). 

It should be noted here that we do not restrict the 
Gaussian operators of this type to be just thermal states. 
In general we would wish to consider density matrices 
that are linear combinations of Gaussian operators, and 
these can exist as hermitian, positive-definite operators 
even when composed of Gaussians that have neither prop- 
erty. 



B. Non- number-conserving Gaussians 

If anomalous products of form bibj are included as well, 
then the most general non-number-conserving Gaussian 
is a type of squeezed state: 



A^^(m,^,I+) = :cxp 



Ul^-^^jbWj] ■■E[^-^^vb^b3 



i>3 

X JJexp 



1 - &^b 



tj"'i-^3 



(2.6) 
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Here ^, are complex antisymmetric M x M matrices. 
For the special case that is a hermitian matrix and 
= then the Gaussian operator is just the density 
operator for a squeezed thermal state. Normalized forms 

Am, for which Tr ^Am^ = 1, will be introduced later. 

As before, we do not restrict the Gaussian operators to 
be just the squeezed thermal states, even though this rep- 
resents a large and physically important class of fermionic 
density operators. By extending the definition to include 
all exponentials of quadratic forms, we can obtain a more 
useful set, which is a complete basis with a positive- 
definite representation for all density operators, as we 
show in Sec. \V\ 



III. ONE AND TWO MODE GAUSSIAN 
OPERATORS 

In this section we give examples of Gaussian operators 
in elementary one and two-mode cases. These are suf- 
ficient to illustrate the basic identities and ideas. More 
general results will be given in Sec. [V] 



A. Single-mode Gaussian operators 

An unnormalised Gaussian operator for a single mode 
has only one possible form: 



exp 



fc=0 



(3.1) 



where the exponential is defined, as indicated, by its se- 
ries representation. Just as in the general case of Eq H2.5|l , 
the anticommutivity of the fermionic operators means 
that only the zeroth and first-order terms in the expan- 
sion contribute to the single mode Gaussian, giving the 
simple result: 



The normalisation of this Gaussian operator is 

Ni = TrAi 
= 2-^1. 



(3.2) 



(3.3) 



Excluding the point /i = 2, we define the new parameter 
n — (1 — /i)/(2 — /i), which allows us to write the Gaussian 
in normalised form as 



Ai{n) = (l-n):exp - (2 + l/(n - 1)) 

= (1 - n)W + nVb 



(l-n)|0)(0|+n|l) (1| 



(3.4) 



If n is real, then Ea H3.4|l shows directly that Ai(n) is the 
density operator corresponding to a mixture of number 
states in the one-mode case. The expectation value of 
the ferniion occupation n = 6^6 is just the variable n. 



1. Completeness and positivity 



As two special cases we obtain the number states: 



|0)(0| ^ Ai(0) 
= Ai(l), 



(3.5) 



which implies that just these two single-mode Gaussians 
form a complete basis set for a number-conserving sub- 
set of Hilbert space. Because super-selection rules pro- 
hibit superpositions of states differing by odd numbers of 
fermions, this is the most general case possible. Further- 
more, we see from Eq ()3.5|) that the most general single- 
mode density matrix can be expanded as a mixture of 
Gaussians with positive coefhcients, since < n < 1: 



p= (l-n)Ai(0)+nAi(l) 



(3.6) 



Additionally in the one-mode case, all physical density 
operators are also Gaussian operators p = Ai (n) , which 
is another proof of positivity and completeness. 

It is clear from these examples that the Gaussian oper- 
ators are overcomplete: a physical density matrix may be 
represented by more than one positive distribution over 
the Gaussians. The most general single-mode Gaussian 
operator, with n complex, provides an even larger over- 
complete basis for physical density matrices, even though 
such Gaussians are not always physical density matrices 
themselves. For example, any uniform-phase distribution 
over Gaussian states gives the zero-particle state: 



|0)(0| 



d0Ai(re**). 



(3.7) 



From the usual hole-particle symmetry of fermion states, 
the single-particle state is similarly obtained from: 



|1>(1| 



d(f>Ai{l-re"''). 



(3.8) 



B. Two-mode number-conserving Gaussian 
operators 

A straightforward extension of the generalized thermal 
Gaussian form to two modes gives 



^2"^ (a*) = : exp 



2J = 1 



1 - ^ Pijh\hj + det/x6|6j&2&i. 



(3.9) 
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where the the elements of the 2x2 matrix fi. The 

last step follows by explicitly expanding the general re- 
sult in Eq 1)2. 5|l . while taking account of the sign changes 
during normal ordering. Again, the series expansion con- 
tains all possible normally ordered nonzero products of 
the blbj pairs. In terms of two-mode number-state pro- 
jectors, the Gaussian operator is 

a(")(/x) = |00)(00| + (l-Ain)|10)(10| 
+ (1-Ai22)|01) (01| 
+ (l-Aiii-/i22+clct/x)|ll)(ll| 
-Aii2|10)(01|-/i2i|01)(10|. (3.10) 



1. Normalisation and Moments 

Following from Eq. (|3.1U|) . the normalisation is 

N2 = 4- 2/iii - 2^22 +det/x. (3.11) 

Defining the matrix n = I — (21 — /x-'") ^, where I is 
the 2x2 identity matrix, we can write the normalised 
two-mode Gaussian as 



A2 (n) = det [I — n] : exp 



-St 21 



b 



= det [I - n] 1 00) (00 1 

+ (nii(l - 7122) + 7112^21) |10) (10| 

+ (ri22(l - nn) + 71127121) |01) (01| 

-h det n 111) (11| + 7121 110) (01| + ni2 |01) (10| . 

(3.12) 

If n is an Hermitian matrix, then the two-mode Gaus- 
sian corresponds to the density matrix of a mixture of 
states of different total number, with coherences ni2 = 
7T.21 between states of the same total number. 

Normally ordered first-order correlations of the Gaus- 
sian density matrices correspond to elements of n: 



Because the matrix n is Hermitian for a Gaussian that 
is a density operator, it can be diagonalised, correspond- 
ing to a change of single-particle basis. In this diago- 
nalised basis, since the coherences are zero, the density 
operator is a mixture of number states, totally charac- 
terised by average occupation numbers. In other words, 
these Gaussian operators correspond to two-mode ther- 
mal states. 



2. Completeness 

We wish to show first that any number-conserving two- 
mode density matrix can be expanded in terms of Gaus- 
sian operators, and second that this can be done with 
positive expansion coefficients. The first result follows 
if we can represent all the number-state projectors be- 
tween states of the same total number using Gaussian 
operators. By inspection of Eq (|3.12|l above, we find 
that: 



1 





1 



"21 



"12 



00) (00 1 


- A2(0) 


10) (10| 




01) (01| 


- A2((^ 


11) (111 


= A2(I) 


10) (01 1 


= 


01) (10| 


= m[ 



ni 
7^2l "2 

ni ni2 

712 



A2( 
A2( 



711 

712 

711 

712 



(3.15) 



Thus, the two-mode Gaussians form a complete operator 
basis for all number-conserving density matrices. 



3. Positivity 



Tr6l6,A2 



(3.13) 



and higher-order correlations reduce to products of first- 
order averages, for example, 



blhblhj^ = 1117122 - 71127121. 



(3.14) 



This kind of factorisation of higher-order correlations 
could be taken as the defining characteristic of a Gaus- 
sian state, and more generally, Gaussian operators, 
rather than the more formal operator definition given 
by Eq. ifT^ . In other words, Gaussian operators have 
both a Gaussian form and generate Gaussian statistics. 
The connection between these two defining features can 
be made more explicit by use of a moment-generating 
function, or characteristic function, which is considered 
in Appendix IbI 



Not only is the two-mode Gaussian a complete repre- 
sentation, but it is also a positive one: any two-mode 
number-conserving density operator can be written as 
a positive distribution over Gaussian operators. To see 
this, note that while the expression for the projection 
operators, Eq H3.15|l . includes terms with negative coeffi- 
cients, the projectors involved are the off-diagonal ones. 
Since density matrices are positive-definite, off-diagonal 
projectors can only occur in combination with diagonal 
projectors. We take this into account in what follows. 

Any two-mode density operator can be expanded into 
number-state projector operators as follows: 



(3.16) 



where n and ft' are vectors of integer occupation numbers: 
n = (711,712), fi' = (n'i,n^). Here pft^n' = if ^ 



5 



J2j n'j , because of total-number conservation. Using the 
relations in Eq. H3.15|l . we can write the density operator 
as 



P 



A2( 



"1 2p(oi).(io) 

712 



) 



+A2( 



(3.17) 



ni 

2P(10),(01) "2 

Since the diagonal coefficients pf^^fi are positive and sum 
to one, the Gaussian operators form the basis of a proba- 
bilistic representation of any two-mode density operator. 

While any two-mode number-conserving density ma- 
trix can be expanded in the form in Eq. H3.17|l . there 
are often more direct representations. For example, 
as Eq. H3.20|l below shows, the entangled state = 
a |10)+/3 |01) can be represented by just one term, rather 
than the four terms that result from Eq. (|3.17|) . 



4-. Correlation and entanglement 

We first note that any uncorrelated product of number 
state mixtures can be represented exactly: 



{(l-ni)|0)(OKni|l)(l|} 
®{(l-n2)|0)(0Kn2|l>(l|} 
ni 

712 



A2( 



(3.18) 



As well as these uncorrelated mixtures, the Gaussian ba- 
sis can also be used to represent a mixture with correla- 
tions between the modes, this time as a sum (with posi- 
tive weights) of two terms: 

^|00)(00|-hB|ll)(ll| = ^A2(0) -hBA2(I). 

(3.19) 

Importantly, superpositions of number states, corre- 
sponding to entangled states, can also be represented, 
subject to total-number conservation. For example, the 
density matrix corresponding to the state \4>) = a |10) -I- 
/3|01) is 

10) (01 = |a|'|10)(10Ka/3*|10)(01| 
+a*/3|01)(10| + |/3|'|01>(01| 



(3.20) 



Thus the number-conserving Gaussian operators even in- 
clude entangled density matrices, making them a power- 
ful tool for representing highly correlated states. 

This property of being able to represent such Bell-like 
entangled states with the hermitian subset of Gaussian 
operators, is different to the case of Gaussian expan- 
sions for bosons. A typical example of this type of non- 
classical representation is the positive P-representation[3 
1^ , which must use a non-hermitian basis to obtain a pos- 
itive distribution that can represent all two-mode densit; 
matrices, such as those that violate a Bell inequality 



isity 

m 



C. Two- mode squeezed Gaussian operators 

Equation (|3.9() does not represent the most general 
two-mode Gaussian form, as the quadratic form does 
not yet include terms such as 6162- Incorporating such 
anomalous products, we can write the most general Gaus- 
sian operator in normalised form as: 



4"'(M,^,r) = :exp 



+ (det/x + ee+)6l65626i. 



(3.21) 



where ^ and ^+ are independent complex numbers. The 
two additional operator terms in the expansion, b\b\ and 

6162, are projectors between states of different total num- 
ber: 



6162 
blh\ 



|00)(11| 
|11)(00| 



(3.22) 



which are the kinds of coherences that appear in the den- 
sity matrices of two-mode squeezed states. 



1. Normalisation and Moments 

The normalisation of the squeezed Gaussian is 

N2 ^ 4-2^11 -2^22+ det/x + ee+- (3.23) 

To incorporate this into a normalised Gaussian, we rede- 
fine the n matrix to be 



n = I 



N2 - ii-^ 



(21 



(3.24) 



and introduce rescaled squeezing parameters 771 ~ 
—^/N2, 771+ = — ^+/A^2- The normalised form is then 



A2(n, 771, 771^ 



det [I — n] + mm'^) x 



: exp 



b 

6tr 



(3.25) 



where the 4x4 matrices / and a are defined to be 



" -1 


0" 








-1 














1 











1. 








"7711-1 


"-21 





777 


?^12 


"22 - 


- 1 


— 771 








— 777 


+ 


1 - 7111 


-"12 


777+ 







-"21 


1 - 7722 



.(3.26) 
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In terms of number-state projectors, the normalised 
Gaussian is a generalisation of the number-conserving 
case, but with the additional non- number-conserving pro- 
jectors: 



over Gaussian operators: 



Mg) = (dot [I - n] + mm+) |00) (00| 

-I- (nii(l - ^22) + rti2"-2i - mm^ 
+ ('^22(1 - '^ii) + «12«21 - mm^ 
+ (detn + TOm+) (11| 
+n2i |10) (01| -Hni2 |01) (10| 
-m|ll) (00|-m+|00)(ll|. 



|10> (101 
|01> (Oil 



(3.27) 



In addition to the normal fluctuations of Eq. I|3.13|) . the 
squeezed Gaussians also contain anomalous fluctuations, 
which are just equal to the new variables m and m+: 



^l&2)~ = ITT" 
A 



(3.28) 



which implies that for A2(n, r7i,m+) to be a density ma- 
trix, m and m+ must be complex-conjugate. The second- 
order correlation generalises to 



6|6i^62)_ = niin22 — ni2n2i + mm^ , (3.29) 
A 



which again corresponds to the decorrelation that occurs 
in a state with Gaussian statistics. 



2. Completeness 

From Eq. (|3.27|) , it follows that the squeezed Gaussians 
provide a complete two-mode fermionic basis, not only for 
the number conserving subspace, but also for all states 
containing superpositions of states whose difference in 
total number is even. To see this, note that the projectors 
between number the |00) and |11) number states can be 
written explicitly in terms of the Gaussian operators as 



TO 111) (00| 
TO+ |00) (11| 



A2(n,-m,0)-A2(n,0,0), 
A2(n,0,-TO+) - A2(n,0,0), 



(3.30) 



for any n. These projectors, together with those of 
Eq. H3.15I) . span the complete Hilbert space of density 
matrices in question. 



3. Positivity 

Just as for the number-conserving case, we can write 
any physical density operator as a positive distribution 



P = 



n n' 
n 



+A2( 



ni 



2p(io),(oi) "2 



, 0, -2p(oo,ii)) 



(3.31) 



Thus the Gaussians form the basis of a probabilistic rep- 
resentation of any physical two-mode density operator. 
The ability to represent physical density matrices with 
the Gaussian basis, either as a single element or by a posi- 
tive distribution over basis elements, is important for cal- 
culating dynamical simulations of quantum systems via 
probabilistic methods. 



4. Entanglement 

Again, these types of Gaussians can directly corre- 
spond to entangled density matrices, without having 
to consider expansions over several elements. But now 
superpositions of different total number can be repre- 
sented. As an example of these additional kinds of physi- 
cal states that the squeezed Gaussians can represent, con- 
sider the entangled non-number-conserving superposition 
state: IV-) = [|00) -f |11)] /%/2. This can be represented 
by a single Gaussian: 

\^){^\ = ^[|00)(00| + |00)(11| + |11)(00| + |11)(11|] 

- A2(il,-i -i). (3.32) 

Hence, the more general type of Gaussian basis element 
considered here is considerably more powerful for repre- 
senting entangled and correlated states than the number- 
conserving basis set. 



IV. MULTIMODE DECOMPOSITION OF 
FERMI SYSTEMS 

So far we have considered only one- and two-mode sys- 
tems, to illustrate the basic physical properties of Gaus- 
sian operators. As we saw, the Gaussian operators could 
easily be written in terms of number-state projectors. 
The power of the Gaussian operators as a basis in its 
own right becomes apparent for multimode systems, for 
which the number-state basis becomes unusable. 

Before defining the general Gaussian basis for a system 
with many degrees of freedom, we define some mathemat- 
ical notation and conventions that will be of subsequent 
use. As before, we define 6 as a column vector of the M 
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annihilation operators, and as a row vector of the cor- 
responding creation operators. We also introduce an ex- 
tended column vector of all 2AI operators: b = {b'^ , b^)'^ , 

with an adjoint row vector defined as fe^ = (6^. b^). Writ- 
ing these out in full, we get: 



where / is the constant matrix 



b = 



DM 



{hi. 



.,b 



M 



(4.1) 



V ill J 

Throughout the paper, we print length-M vectors and 
M X M matrices in bold type, and index them where 
necessary with Latin indices: j = 1,...,M. Length-2M 
vectors we denote with an underline and 2M x 2M matri- 
ces we denote with a double underline. These extended 
vectors and matrices are indexed where necessary with 
Greek indices: /i — 1, ...,2M. Note that an object such 

as & 6 is a 2M x 2M matrix: 



bb 



bbt bb^ 

bt^bt bt^b^ 



(4.2) 



whereas b b = J2k + ^) ^ scalar. More 

general kinds of vectors are denoted with an arrow nota- 
tion: ~jj . 

For products of operators, we make use of normal and 
antinornial ordering concepts. Normal ordering, denoted 
by :•••:, is defined as in the bosonic case, with all anni- 
hilation operators to the right of the creation operators, 
except that each pairwise reordering involved induces a 
sign change, e. g. : : = — • The sign changes 
are necessary so that the anticommuting natures of the 
Fermi operators can be accommodated without ambigu- 
ity. We define aniii normal ordering similarly, and de- 
note it via curly braces: {b'jbi} = • More gener- 
ally, we can define nested orderings, in which the outer 
ordering does not reorder the inner one. For example, 
{: 05j : bi} — —bib'j : O : , where O is some operator. 

For example, the different orderings of pairs are, in 
block matrix form. 



bf,bl : = - : blb^ : = 



-bjbi bibj 
bib] 

hb^j bib J 

W] -bM 



(4.3) 



Note that this convention means that the relation be- 
tween the two orderings is 



:bb : = 



L+{bV}, 



(4.4) 



/ 
/ 



(4.5) 



in which and I are the M x M zero and identity ma- 
trices, respectively. 

When ordering products that contain a Gaussian oper- 
ator A (and later the density operator), we do not change 
the ordering of A itself; the other operators are merely 
reordered around it. Thus |: A6| : 6j| = — 6j6|A. The 

different possible quadratic orderings containing a Gaus- 
sian operator can thus be written in matrix form as 



bbA: = 



{bVA} 

{b-.n-.} 

{■■ A6 : V} 



6Abt SgTA 



AbtTfot 



bb^A 

-(stTAfety 

AbSt 
St^Agt 



bAb^ 
bbtA 



bAb^ 
Abbt 

T 



(4.6) 



V. NORMALIZED GAUSSIAN OPERATORS 

A. Definition 

We define a normalized Gaussian operator A to be 
the most general Gaussian form of fermionic annihila- 
tion and creation operators, with zero displacement and 
unit trace. Using the extended-vector notation, we first 
write the normally-ordered Gaussian in an unnormalized 
form as 



A("H/^,$,I+) = :exp V {L^g-'/2)b 



(5.1) 



where we have introduced a new extended 2AI x 2M co- 
variance matrix a defined in terms of /x, ^, ^+ so that: 



o-~i = 21 



(5.2) 



The introduction of the generalized covariance g_ allows 
the matrix to be written in a normalised form, using the 
results of Appendix IBl together with an explicit complex 
amplitude factor Q: 
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a(a: 



exp 



F (l-g~V2)& 



(5.3) 



The normalisation is one obvious difference with the con- 
ventional complex- number or bosonic Gaussian forms Q. 
Chosen to ensure that Tr A — ft , it contains the Pfaf- 
fian of an antisymmetric form cta of the covariance. The 
choice of anti-symmetrisation is given in Appendix ^ 
other choices will lead, in general, to additional sign fac- 
tors. Now the square of the Pfafhan of an antisymmetric 
matrix is equal to its determinant, and the determinant 
of (TA differs from that of a; by a constant sign (see Ap- 

pendixEl). Thus |Pf 



det I 



and as we shall 



see, the relative phase between the two does not appear 
in later identities. The additional variable il plays the 
role of a weighting factor in the expansion that allows 
us to represent unnormalised density o per ators and to 
introduce 'stochastic gauges' in the drift 21,'25|. 

The covariance has a type of generalized Hermitian 
antisymmetry, which can be written as = — with 
the definition that: 



a 


b 


+ 


■ d 


c 


c 


d 




b 


a 



(5.4) 



It is this generalised antisymmetry that allows the covari- 
ance to be transformed into an explicitly antisymmetric 
matrix. The covariance can also be broken down into the 
physically significant M x M submatrices n , m and m+ : 



m ' 



— I m 



(5.5) 



Here n is a complex matrix, which corresponds to the nor- 
mal Green's function in many-body terminology, while m 
and m+ are two independent antisymmetric complex ma- 
trices that correspond to anomalous Green's functions, as 
we will show in the next section. 

Thus the total parametrization of a general Gaussian 
operator is 



A — {n, n, m, m+) , 



(5.6) 



consisting of 1 -I- p = 1-1- M{2M — 1) parameters in all. 
However, for many systems, only a subset of all Gaus- 
sian operators is required for a complete representation 
of the density operator. One important subset is the set 
of generalised thermal states, for which m = m+ = 0. 
In this case, from Appendix IBl and using the result that 
[21— /i.] = [I— n]"-'^, the normalization factor is det [I — n]. 
The normalized thermal Gaussians therefore can be writ- 
ten: 



a(a; 



n det [I - n] : exp -b^ (21 + (n'^ - I 



In order to use the Gaussian operator basis, we need to 
make use of a number of basic results. The proof of many 
of these can be established, as we show in Appendix IbI 
with Fermi coherent states and Grassmann algebra, the 
basics of which are given in Appendix ^ However the 
final results do not contain any Grassmann variables. 



B. Trace Properties 



(5.7) 



Some basic traces are 





Tr 


A 






Tr 


bA 


= , 


Tr 




A : 


= Qa 



(5.8) 



The first of these is the normalisation, proved as theorem 
1 in Appendix ^ That the second is zero follows from 
the fact that the Gaussians are constructed from pairs of 
ladder operators and thus cannot correspond to a super- 
position of states whose total fermion numbers differ by 
one. The same result holds for the trace with any odd 
product. The third trace, proved as theorem 2 in Ap- 
pendix^ allows us to calculate first-order moments. In 
terms of the M x M submatrices, these become: 



Tr bjb.A 



= fir 



Tr 
Tr 



bibj A 



(5.9) 



These results imply that for A itself to correspond to a 
physical density matrix, n must be a Hermitian matrix, 
t;, .\ — /;^t;,.\ ^j-^^ jjjyg^ ]-,g |;[jg Hermitian 



since (blbj 



b]b, 



conjugate matrix of m, since {bibjj — {b^b-j. Physi- 
cally, n gives the number, or normal, correlations, and m 
and m"'" give the squeezing, or anomalous, correlations. 

Another restriction on A being a physical density ma- 
trix that follows from Eq. (|5.9II is that the eigenvalues 
of the matrix n lie in the interval [0, 1], because of the 
Pauli exclusion principle for fermions. Furthermore, the 
variance in the number correlations is 

var{(6t6,)} = (bib,) - (Wj)) ^ (5-10) 

which implies that if all the eigenvalues of n are or 
1, then A is a number state, as the variance in number 
vanishes. If the eigenvalues are not limited to or 1, 
then the A corresponds to a mixture of number states in 
the eigenbasis and is thus a (possibly squeezed) thermal 
state, characterised by average occupation numbers nj = 
eigj(n), and squeezing matrix m. 
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A general Gaussian operator will not necessarily ful- 
fill the Hermiticity condition and thus will not neces- 
sarily correspond to a physical state. However the set 
of operators that do correspond to physical states is an 
important subclass, because the expansion allows these 
states to be represented with exact precision. The inclu- 
sion of nonphysical operators in the expansion, on the 
other hand, makes the Gaussian basis an (over)complete 
basis in which to expand a physical density operator of 
an arbitrary state (This is proved below for the general 
case). The overcompleteness of the Gaussian operators 
as a basis set has important implications for representing 
arbitrary states with a positive distribution function, a 
fact that we discuss in detail elsewhere [13. 



C. Completeness and Positivity 

We next wish to show that the previous results on com- 
pleteness and positivity obtained for one and two mode 
density matrix representations can be generalized to the 
multi-mode case. That is, we will prove that: 

• For any physical density matrix p, a positive set of 
coefficients Pj exists such that 

3 

This central result does not rely on utilising the com- 
plex amplitudes fi, which are part of the most general 
Gaussian operator. If these were used, then positivity 
of the coefficients would be trivial, since these additional 
amplitudes could be used to absorb any phase or sign fac- 
tors arising in the expansion. Instead, we wish to prove 
a much stronger result, that a positive expansion exists 
without any additional amplitude factors. This result is 
analogous to a similar result known for the positivc-P 
bosonic representationj3,l3. 

From the number state basis of Eq (|2.2(l . the full set 
of possible fermionic many-body number states is the set 
{\n)} where n is varied over all 2^^ possible permuta- 
tions. This defines a complete basis of 2^*^ number-state 
projectors for the set of all fermionic density operators. 
While not all the number-state projectors are Hermitian, 
it is no restriction to use this larger set of operators as a 
basis for the density matrices p. 

Next, we expand: 

P = I") {^\p\rn) {rn\ 

n rh 

= X! XI 1^) 

n rh 
n m 

The positive definiteness of the density operator means in 
particular that diagonal density matrix elements are real 



and positive: pftft > 0. It is sufficient for completeness 
to prove that any elementary fermionic operator of form 
Piirn = Prim \n) (™| corresponds to a normalised Gaussian 
A( A ), apart from a positive scaling factor. The demon- 
stration proceeds by constructing limiting cases of Gaus- 
sians that correspond to each of the elementary compo- 
nents of the density matrix. As demonstrated in the one 
and two mode cases, such expansions are not unique, and 
generally one can obtain other more compact representa- 
tions by combining diagonal and off-diagonal elements. 

To prove the elementary result, we proceed in three 
steps: 

1. Diagonal operators 

First, the generic diagonal operator in the number ba- 
sis is: 

Pnn = Prm \n) (n| 
i 

with a total occupation number 

iV^^n, . (5.14) 

i 

Each diagonal multi-mode projector \n) {n\ is simply an 
outer-product of single-mode density matrices. From the 
single-mode example in Section Fill Al each single-mode 
density matrix corresponds to a Gaussian as in Ea H3.6|l . 
Thus, we see that a diagonal projector is exactly equal to 
a normalized Gaussian (we suppress the trivial arguments 
for simplicity): 

\n) (n| = n. Al K] 

= Am [n] , (5.15) 

where the matrix n is defined as — riiSij . 

Because the pfifi are real and positive, Eq. H5.13|l shows 
that a positive Gaussian expansion exists for all diagonal 
density matrices: 

Pnn = Pnn^M [n] . (5.16) 

In summary, a diagonal multi-mode projector \n) {n\ is 
simply an outer-product of single-mode density matrices, 
and hence corresponds exactly to a multi-mode normal- 
ized Gaussian. 



2. Generalized thermal operators 

Next, consider off-diagonal projectors in the number- 
state expansion that conserve total number, i.e. 

Pnni = Pnrh \n) (m| , (5.17) 
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for which 



N 



(5.18) 



We show that any such component can be written as the 
limiting form of a number-conserving Gaussian, up to a 
positive scahng factor. 
Let 



n- = (5y minjni, mj, 



(5.19) 



and define an off-diagonal Gaussian in terms of the diag- 
onal normalized Gaussian: 



A[n'] 



(5.20) 



Now for every distinct mode i with Sui = rii ~ 7i[ = 1 
one can define a corresponding distinct index with 
5m j = rrij — m ^ = 1. It follows from the minimum con- 
dition, Ea l5.19| ) that ^ i, and from the number- 
conservation equation, Eq H5.18|l . we have 



(5.21) 



Similarly, for every distinct pairs of indices with 
Sm = Suii = 1, it follows that ^ i, since other- 

wise Sui = 0. The mapping therefore generates distinct 
pairs so that ^ j(i') <^=^ i ^ i' . Next, we note that 
this mapping is not a permutation of the set of modes i 
with 5ni = 0, since if it were the condition that ^ i 
would be violated for some i' . Similarly, the mapping is 
not a permutation of any subset of these modes. This 
means that the only non-vanishing terms in the normal- 
ization factor det[2I — fi] are the diagonal terms, which 
are already normalized. 

Proceed by defining the resulting set of SM < M pairs 
as a = , and let: 



^Sii' 5. 



{i',f}Ga 



where ii = pnm, so that: 

A[/x(£),n'] = -.Ull-fMibrbj 



A[n']: 



n 



1 



(5.22) 

Then consider the limit of e — > 0, so that to leading order, 



Pnri 



SM 



n p 



lim £"^A[/x(£),n'] 

£— >0 



(5.23) 



Again, we see that a positive expansion parameter is 
obtained. 



3. Squeezed operators 

Finally, we consider the remaining elements of the den- 
sity operator expansion, i.e. those squeezed projectors 
'Pfim for which 



25' 



(5.24) 



where S is an integer denoting the change in the number 
of fermion pairs. We suppose that 5 > 0, since the case 
of S' < is obtained trivially by hermitian conjugation. 
Let rii be obtained from Ui by removing 25* occupied 
sites, labeled as successive pairs z, j belonging to a set tr, 
so that N = X^i'^j ~ "The occupation numbers 

ni,mi now define a generalised thermal density matrix 
component as previously, and hence equate to a limiting 
case of a Gaussian operator from Eq (|5.23() above. 

Now define a squeezed off-diagonal Gaussian in terms 
of the thermal case, which we have already proved has a 
positive representation: 



A(^)(m(£), ^(e), n') = n [l - ^^Mbl\ A [^(e), n'] 



where: 



E 



— £ 



(5.25) 



(5.26) 



Then consider the limit of £ —> as before, so that to 
leading order. 



Pnrn 



n M 



= lini£^+^'*'^A(^)[/x(£),|(£),n'] (5.27) 

A positive expansion parameter is obtained here as 
well, thus completing the proof. 



D. Differential Properties 

In order to use the Gaussian basis in a time-evolution 
problem, we need to be able to map the evolution of the 
density operator onto an evolution of the expansion coef- 
ficients Pj. To achieve this, one must be able to write the 
action of ladder operators on a Gaussian basis element 
in differential form. 

We can differentiate the Gaussian operators with re- 
spect to their parameters to get 

— A = —A, 



da — 



(5.28) 
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where the matrix derivative is defined as 

d 




(5.29) 



i. e. involving a transpose. These expressions for the 
derivative can be inverted to obtain the important iden- 
tities: 



A 



bb A : 



T dA 
a A — a — — a . 
— — da — 



(5.30) 



and thus we can write the normally ordered action of any 
pair of operators on a Gaussian as a first-order deriva- 
tive. As theorems 5 and 6 in Appendix iBlshow. there are 
analogous identities for antinormally ordered and mixed 
pairs: 



7 ^ r\9A 



{b-.n-.} 

{bVA} = {g-l)A-{g-l)^{g~l). 



(5.31) 

For the subset of Gaussian operators that correspond 
to (generalised) thermal states, i. e. m+ = m = 0, we 
obtain a simpler set of differential identities: 





= nA 


,dA 
HI-n)^n. 






^ nA 


^n-il-n) 






- (I- 


- n)A+ {I - n) 




(bApy 


= nA 


dA 
- n—n. 
an 


(5.32) 



The action of four ladder operators on a Gaussian op- 
erator can be obtained by applying the previous identi- 
ties twice. Thus in a Gaussian expansion, any two-body 
interaction term can be written as a second-order differ- 
ential operator. As we show in detail elsewhere(15j. this 
central result allows the evolving density operator to be 
mapped to a Fokker-Planck equation for the expansion 
coefficients, thus enabling a Monte-Carlo sampling of the 
many-body quantum state. 



special cases, the set of Gaussian operators include the 
density operators for thermal states and squeezed states, 
and thus the physics of the noninteracting Fermi gas and 
the BCS state is incorporated into the basis itself. Fur- 
thermore, the basis also includes more general operators, 
which ensure an overcompleteness that makes it possible 
to express any physical density operator as a probabilis- 
tic distribution over the Gaussian operators, without the 
need to use Grassmann algebra. 

We have calculated the normalisation and moments, 
and proved completeness and positivity for the most 
general basis, and also for specific subsets, such as the 
number-conserving thermal basis. These results mean 
that the phase-space representation defined by the Gaus- 
sian operators can be used for first-principles simulations 
of many-body fermionic systems. The mapping from the 
quantum operator to the probabilistic c-number descrip- 
tion is enabled by the Gaussian differential identities that 
have been derived here. The application of these identi- 
ties will be dealt with elsewhere, when we consider the 
phase-space representation in more detail. 
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Appendix A: GRASSMANN ALGEBRA 

This appendix introduces the basic concepts of non- 
commuting algebra and lists some results pertaining to 
Grassmann calculus and Fermi coherent states. These 
results are used in Appendix ^ to establish important 
properties of the Gaussian operators. Proofs and fur- 
ther discussion of these results can be found in the 
literature_16, 17, 18]. 

Let q; be a vector of M anticommuting (Grassmann) 
numbers, i.e. 







(Al) 



Since = 0, any function of Grassmann numbers can 
be at most linear in any one of its arguments. Thus, for 
example, the single mode exponential is 



exp(aj) = 1 



(A2) 



VI. CONCLUSION 

In summary, we have introduced here a generalised 
Gaussian operator basis, as a means of defining a phase- 
space representation for correlated fermionic states. As 



and a multimode exponential, e. g. exp (^ aj) will be an 
ordered product of such single-mode exponentials. 

The Grassmann numbers anticommute with all Fermi 
annihilation and creation operators, but commute with 
c-numbers and bosonic operators. 
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Grassmann Calculus 



Grassmann coherent states 



Differentiation of a single Grassmann variable is de- 
fined as 



A - A 



(A3) 



with the derivative of products obtained by the Grass- 
mann chain rule: 

i. e. the derivative operator also anticommutes. 

Grassmann integration is defined to be the same as 
differentiation, but written in a different way: 



/ 



daj = 



1 . 



(A5) 



Multivariate integrals are ordered sequences of single- 
variable integrations, which can be written without am- 
biguity if the integration measure (as for a derivative) 
is also taken to be an anticommuting number. For an 
integral over all variables in a vector, we define the in- 
tegration measure to be ordered in increasing numerical 
order: da = dai-.-da^ ■ 

Note that integrating a total derivative gives zero: 



. 



(A6) 



This fact, coupled with the product rule, gives a result 
for partial integration: 



a/ _ \-!da,f^ , /even 
/ 1 Jda.f^ , /odd • ^^'^ 



One very useful result is the multivariate Gaussian in- 
tegral: 



J exp{-a^Aa/2)da = Pf[A] 



(A8) 



for an antisymmetric matrix A of complex numbers. The 
PfaSian is related to the determinant (Pf [A])^ = det A , 
and thus, apart from a sign change, has many of the 
same properties. For example, Pf [A^] = (— l)''^Pf[A] 
and |Pf [A-i] | = |l/Pf[A]|. Another useful type of 
Gaussian integral is: 



/ 



M 



exp (-/3^Ba) J| {d^jdaj) 



det [B] , (A9) 



where B is a square matrix of complex numbers and /3 
and a are two independent Grassmann vectors. This 
second integral is in fact a special case of the first, except 
with 2M total Grassmann variables. 



For each Grassmann number aj we can associate an- 
other Grassmann number, denoted aj , to play the role 
of a complex conjugate. This is formally regarded as an 
independent variable for calculus purposes. The conju- 
gates aj anticommute with all other Grassmann vari- 
ables. By use of such a complex Grassmann algebra, we 
can define a fermionic coherent state, which is formally 
an eigenstate of the annihilation operator: 

\aj) = (l + a,a,)-M|0) + |l)aj) 
= exp (b]aj - ^ajaj^ |0) , 



(AlO) 



{aj\ = {l + a,a,) M(0|+S,(1|) 
= exp (ajbj - ^ajttjj (0| . 



(All) 



Like the bosonic coherent state, the fermionic coherent 
state can be written as (Grassmann) displacement from 
the vacuum: 



exp 



(S]a, -aA) |0) 



(A12) 



Multimode coherent states are products of the single- 
mode states: 



M 



i«) = n 1"^) 



= exp 



M 



Y^b]aj - ajbj 



|0) 



The inner product of two states is: 

{ai\ aj) = exp (a,Q!j — aiai/2 — ajaj/2) , (A13) 

and thus as two special cases: 

{ai\ai) = 1 
{-ai\ai) = exp(-2aiai) . (A14) 

The usefulness of the coherent states lies in the fact 
that they form a complete set: 



dajdaj \aj}{aj\ = |0)(0| + |1) .(1| = h , (A15) 
or, for the multimode case, 

J da\a){a\ = Im , (A16) 

where we have defined a 2M-variate integration measure 

as 



M 



da = JJ {daj daj) . 



(A17) 
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Finally, we can express the trace of an arbitrary operator 
O as 



Tr 



O 



= da (—ot\ O la) 



(A18) 



Appendix B: GAUSSIAN FERMION OPERATORS 

We prove some useful results concerning the most 
general multi-mode Gaussian operator constructed from 
fermionic ladder operators. The proofs make use of the 
properties of Grassmann coherent states and anticom- 
muting algebra, which are summarized in Appendix 1X1 
However the final results do not contain any Grassmann 
variables. The results establish the trace and differential 
properties of the Gaussian operators that are discussed 
in Sec. El 



1. General Gaussian Operator 

In this appendix, we use the vector and ordering nota- 
tions introduced in Sec. lIVI We find it convenient to use 
an unnormalised Gaussian form of Fermi operators: 

A(")(^) = : exp :, 

/ 1^ ^ \ 

= n 1- 2^^.Mm-^^. (Bl) 
where /i is a 2M x 2M matrix of parameters, given by 



In terms of the covariance matrix a, the unnormalised 
Gaussian is: 



A(")(o:) 



exp 



(B2) 



where the relation between the two parametrizations is 
gi — + 2J]~^ and where the diagonal matrix / is as 

defined in Eq. H4.5|l . 

Because of the anticommuting property of fermions, 
both a and /x possess a generalized antisymmetry: a = 

-fJ,'^ , or more in block matrix form. 



a b 
c d 



= -X 



a b 
c d 

b^ 
a^ 



X 



(B3) 



where the constant matrix X is defined as 



X 



I 

1 



(B4) 



When applied to the left of a matrix, X swaps the upper 
and lower halves; when applied to the right, it swaps 
the left and right halves. This structure means that 
the matrices a and fj, can be transformed into explic- 
itly antisymmetric forms by certain permutations of rows 
and columns. For example, it follows immediately from 
Eq. HB3|I that interchanging left and right halves, or up- 
per and lower halves, will generate an antisymmetric ma- 
trix. Alternatively, inserting each row in the lower half 
after the corresponding row in the upper half, and insert- 
ing each column in the right half before the corresponding 
row in the left half generates the antisymmetric form 



a b 
c d 



hi 
dii 



ail 

Cll 



bMl flA/i 
dMi cmi 



i>iM aiM 

dlM CiM 

bMM aMM 

dMM cmm 



.(B5) 



With the covariance matrix antisymmetrized in this way, 
the Gaussian operator becomes 



exp 



£(l^-^^V2)L] : , (B6) 



where the vector of operators is now 

(hi \ 



(B7) 



The generalised antisymmetry of a and fi has impli- 
cations for matrix derivatives. Because each element 
of the matrix appears twice, we have da^iv / dag^ = 
Sfj.e^i^4> — X^^X^g, where 5^^ is the Kronecker delta func- 
tion. The extra terms here appear in the derivative of an 
inverse: 



^ Xi4>^0ccr ^J- , (B8) 



and they also give an additional factor of two in the 
derivative of a determinant, thus: 



det (7 (7 ^ , 



— 1 / det a 
da 



where we define the matrix derivative as 

d 




dOyn 



(B9) 



(BIO) 
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i. e. involving a transpose. The result IjBSp allows us to 
relate the derivatives with respect to g_ and /i : 



df 



df 



da ^ — da — 



Normalisation 



(Bll) 



Theorem 1: The trace of the unnormalised Gaussian 
operator is equal to the Pfafhan of the inverse of the 
antisymmetrized covariance, i. e. 



N =Tt 



A(")(<t) 



Pf 



(B12) 



where Pf stands for Pfafhan. 



Proof: Because the Gaussian is in normally ordered 
form, it is straightforward to evaluate the trace with mul- 
timode Grassmann coherent states \o^)^ using the Grass- 
mann trace result of Eq. (| A 18(1 : 



Tr 



M 



3 

» M 

/ cxp [— q_^^ W [—dxxjdaj] , 

(B13) 



where we have made use of the fact that, from the Grass- 
mann inner product result of Eq. I|A13|) . (— a] ex) = 
exp(— 2aa) . We have also changed variables: olj 
—olj , and introduced 2Af- vectors of Grassmann variables 
a — (a, {pt)^) and = (a, a^) . Changing to the anti- 
symmetric form of the covariance a^, and swapping the 
order of the pairs in the integration measure, we obtain 



Tr 



A(") 



M 



exp 



T 



(B14) 



where the reordered Grassmann vector is = 
(ai, ai, aM ,oIm)'^ ■ Noting that the arrangement of el- 
ements in matches the order of integration, we can 
apply the Gaussian integral result [Ea. lA8|l ]: 



Tr 



A(^ 



Pf 



(B15) 



QED. 



Corollaries: 

Now the square of the Pfafhan of a matrix is equal to its 
determinant. Thus, to within a ± sign, the normalisation 



is determined by the determinant of the covariance: 



1 



1 



det 



det det [i^ 

(B16) 



and thus we may write the Gaussian operator in nor- 
malised form as 



A = ±y'det[i£] :exp f {L-g^^f^)! ■ , (B17) 
where the choice of plus/minus sign is determined by 



det 





Pf 


'gA\ 









This extra sign, or phase, which 

does not appear in the normalizations of the familiar 
complex- number or bosonic Gaussians, fortunately does 
not appear in any of the identities needed to make use 
of the Gaussian operators as a basis for a phase-space 
representation. 

A specific case where the determinant appears is for 
the generalised thermal Gaussian without squeezing pa- 
rameters, so that m = m+ — 0. In this case, the nor- 
malization follows directly from the second Grassmann 
Gaussian integral identity, Eq. (jA9|) . Following the same 
procedure as before, we find that: 



Tr [A(")j = det [21 - /x] 



First-order moments 



(BIS) 



Theorem 2: The fermionic Gaussian operator is com- 
pletely characterised by its first-order moments. In par- 
ticular, the covariance matrix corresponds to the first- 
order moments in normally ordered trace form, i. e. 



Tr 



bb A : 



(B19) 



where bb is a matrix multiplication of two vectors, re- 
suhing in the 2M x 2 A/ matrix of Eq. (|0) . 
Proof: We proceed as in the proof of Theorem 1, by tak- 
ing the trace of the unnormalised form using Grassmann 
coherent states and then changing variables: olj — )■ —aj : 



Tr 



66 A(" 



M 



-cx\ ; 66^A(") : \a.) \[ [dajduj] 

3 

» M 

= / aa^ exp [— a^ff "'^ a/2] JJ^ [da.jdxxj\ . 

(B20) 

Next we put the integral into the form Eq. I|A8|I : 
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Tr 



J — (a^"^ a"^)''" exp [—a} a_ ^ a/2] [dajdaj] 

J 

-j—^ / exp [-a'gr^a/2] J| [da^da^] , (B21) 



4. Higher-order moments 

One can calculate higher-order moments along similar 
lines, i. e. by expanding the trace as a Grassmann integral 
then converting this into a higher-order derivative of a 
determinant. The results of this procedure in the general 
case can be written in terms of a moment generating 
function. 



where in taking the derivative with respect to a_~^ we 
have taken account of the fact that each element of it 
appears twice, owing to its generalized antisymmetry. 
Evaluating the Grassmann integral, and employing the 
determinant result , we get 



Tr 



: 65 A(") 



= ±Jdet\i^-^g . (B22) 



Finally, dividing through by the normalisation in 
Eq. I|B12|I gives the normalised result . QED. 

Corollaries: 

We can put Eq. (jHTl) into a more familiar form by 
using the cyclic property of trace: 



Tr 



bb^A 



= Tr 



= Tr 



- (bti^Ab^) Abb^ 
bt^btA bt^^Ab^ 

bbt bb^ 



bt^^bt (bbt 



A 



n I m 

m+ I n 



(B23) 



where we have defined the matrix n for the number, or 
normal, moments, and the matrices m and m"*" for the 
squeezing, or anomalous, moments, as follows: 



n = /bt^b^)„ , 



m 



bb 



/a 

m+ = /bt^bt 



where ( O ) _ = Tr 



OA 



(B24) 



Thus we can write covariance 



matrix in terms of the moments as 

I m 



a 



n 



m ' 



I-n 



(B25) 



Theorem 3: Any even moment of a Gaussian opera- 
tor, in normally ordered trace form, can be calculated by 
means of the moment generating function 



M(x) = V^et[l2_g^ 



(B27) 



as follows 
Tr 



-M{t) 



T=0 

(B28) 



Proof: We start by writing the normally ordered oper- 
ator product as a derivative of a normally ordered Gaus- 
sian: 



d 



d 



: exp 



r=0 



(B29) 



where r is a 2M x 2AI matrix of complex numbers with 
the generalised antisymmetry r = — t+ . Using this result 
in Eq. (|B28() shows that the moment generating function 
M(t) must satisfy 



M{t) = Tr : exp [b t6/2) A 



(B30) 



which can be evaluated by writing the trace as a Grass- 
mann integral 

„ M 

^fe) =7! /exp [-(^ {gr^ - z) a/2] Y[ [dajdaj] . 



Using the Gaussian integral result 



we get 



M{t) = Pf 



-La 



/Pf 



(B31) 



(B32) 



Because of the relationship between Pfaffians and deter- 
minants, we can rewrite this as 



or inverting the relationship. 



: 66 :)„ = I-Ial. 

A 



(B26) 



M{t) 



det 



/,/det 



det [l^ - gz\ , 



(B33) 
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where the sign of the square root is chosen to give a 

positive result when r = 0. QED 

Examples: 

To calculate the derivatives of the moment generating 
function, one makes use of the results for the derivative 
of a determinant Eq. (|B9p and of an inverse Eq. (|B8|I . A 
general second-order moment is of the form 



Tr 



' til p,3 1-1.2 

(B34) 



Thus the normally ordered number-number correlations 



Tr 



: n,;Tij : A 



rriijirnt, (B35) 



i. e. containing the three terms expected for a state with 
Gaussian correlations. 

Similarly, a third-order moment is of the form 



Tr 



^ ll -l II.A ^tl.^ II.R ^ 



'/JlM2'^M3M4^M5A'6 ^ ^MlM4'^M3/J6'^M5M2 ~ "'mi M6 ^M3 ^"2 "'ps 

+ C^M3M6 (^^)piP5 (^^)m4A'2 ^ '^A'5M2 (c^^)piM3 (^^)p6P4 
+ '^.5.4(^^),,,3(^^).e.2 ■ (B36) 

Thus the normally ordered triple correlations are: 



Tr 



UiUjUk : A = n.iTijjnkk - nu {rijknkj + rrijk'm'^f, 

+ nijUjknki - Ujj (uikUki + mikinff.) 

+ njiUkjUik - rikk {nijUji + rriijmf^) 

+ ntjiTiikm^f^ + njimjkmlf. 

+ njkm.jim'l^ + nkjnikimt 

+ nkinikjmfj + nikmkjm^^, (B37) 

again as expected for a state with Gaussian correlations. 



5. Normally ordered products 



Theorem 4- A normally ordered product of a pair of 
ladder operators and a Gaussian is equivalent to a first- 
order differential operator on the Gaussian, as follows: 



bb A : = crA - a— a 
— —oa — 



(B38) 



Proof: The proof can be established easily without 
Grassmann algebra. We first take the derivative of the 
unnormalised Gaussian operator: 



d 



A(«)(m) =- :blA(-\^,)K: . (B39) 



We write this as a derivative with respect to the covari- 
ance matrix, using Eq. HB11|I . and swap the pair of oper- 
ators, to give, 

a T-^q_ ^- :&& A(")(o:): . (B40) 

— da — — 

Next, we take the derivative of the normalisation, using 
Eq. |B9^ : 

;iA' = >)(d..[.d)- 



(B41) 



Combining both of these results, we find that the deriva- 
tive of the normalised Gaussian is 

^A = -iV-2A('0^?^+iv-ii.A(") 
da da da 



= a-^A-a-^ -.bVl : a-\ 



(B42) 



whose inverse is the required result. QED. 

This result can also be proved by use of Grass- 
mann coherent-state expansions, in similar manner to the 
proofs below for the products of different ordering. 



6. Mixed products 



Theorem 5: A product of mixed order of a pair of ladder 
operators and a Gaussian is equivalent to a first-order 
differential operator on the Gaussian, as follows: 



{b-.n-.} 



dA 

-a A +(a-l) —a . (B43) 
— ^— — ' oa — 



Proof: We first make use of the Fermi coherent-state 
completeness identity to replace all ladder operators by 
Grassmann integrals over coherent projection operators: 

{S:S^A:} = 

d7d^dade|7)(7| : |/3> (/3||^A |a) (a| :} |e>(e| 

:— J d'yd(3dade\j){e\exp 
/3 



-(/3,a) (m + I) ( ^ )/2 



_ ] (/3, a) exp 



7/3 + ae — aa ~ (3(3 — -^11 ^ 

(B44) 
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where we have used the result that the inner product of 
two coherent states is, from Eq. (|A13() 

{(3\a.) = ex:p\poL-^l3/2-aoL/2\ . (B45) 

Next, we employ integration by parts [Eq. l|A7|l ] to 

'^^ 

OL 

form: 



replace ( ^ ) by variables that appear in the Gaussian 



exp 



= exp 



^ ) exp {-aa. - (30) 



{13, a) 



-(/3,a) + ( ^ )/2 



exp (— — /3/3) 



(/3,a) 



exp [—oca — (3/3^ 



-dr. 



exp 



+ ( ^ )/2 



cy. 



(/3, 



7. Antinormal products 



Theorem 6: An antinormally ordered product of a pair 
of ladder operators and a Gaussian is equivalent to a first- 
order differential operator on the Gaussian, as follows: 



{fci^A} = ig-L)^- {g-L)^{g-L) (B49) 



Proof: The proof initially proceeds in the same man- 
ner as for products of mixed ordering. We first insert 
the coherent state identity to convert the action of the 
operators into integrals over coherent-state projectors: 



+ ij J d2dpdadeh){'r\ {bP |/3)(/3| A \e){e\ 



= exp (^—aa — (30) 

xexp - (Aa) + (^^ j /2 

(B46) 

We can now express the result as a derivative of the 
unnormalised Gaussian operator with respect to ^ : 



N 



djdf3dade I7) (e| exp 



-(/3,«) (m + I) ( ^ )/2 



X I ^ 1 (q;,/3)cxp 



7/3 + ae — acx — (3(3 — -^1 ^ -^^^ 

(B50j 



|| : S^A :} = ^ y" djdddade \j) { 



This time, however, we integrate by parts twice: 



X exp 



aa — (3(3 — t^I — ^ '^^ 



exp 



exp 



(A«)(^ + l)(^)/2 



1 

'TV 



^ ) (a, /3) exp {-OLOL - (3(3) 



(B47) 



: exp 



(/3,a)(^ + l)(^)/2 



Finally, we can change variables to g_= [/i + 2/] ^ and 
use the result for the derivative of the normalisation: 



^ ) {da, ~c^) exp {-aoi - (}(}) 



{b-.n-.} 



1 

N 



-L{g-' - r) g^g- LA^''^ 



■exp{-aa- (3(3) ^ 



X exp 



-{(3, a) (fi + l 



/2 



= exp (— ao; — /3/3) 



-crA 



. dA 



(^ + l)lexp ~{(3,a){g + lj(^]/2 



(B51) 



(B48) 



QED. 



which is now in a form that we can again express as a 
derivative of the unnormalised Gaussian operator with 
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respect to /i : 



{M^a} = ^ j d2dpdade\-f){e 



X exp 



cxcy. — ^f^f^ ^ 77 ^ + 7/3 + 



dfi 



X U + I I 



1 



-/A(") +l(n + I 



exp 



djji 



Finally, we change variables to ct = [At + 27] -'^and use the 
result for the derivative of the normalisation: 



dk 



(g-l)A-(g-l)-(g-l) 



(l-g- 
(Bf 



(B52) 



QED. 
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